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Abstract 

We present a numerical study of the dynamics of a non-ideal fluid subject to 
a density-dependent pseudo-potential characterized by a hierarchy of nested 
attractive and repulsive interactions. It is shown that above a critical thresh- 
old of the interaction strength, the competition between stable and unstable 
regions results in a short-ranged disordered fluid pattern with sharp density 
contrasts. These disordered configurations contrast with phase-separation 
scenarios typically observed in binary fluids. The present results indicate 
that frustration can be modelled within the framework of a suitable one-body 
effective Boltzmann equation. The lattice implementation of such an effective 
Boltzmann equation may be seen as a preliminary step towards the devel- 
opment of complementary /alternative approaches to truly atomistic methods 
for the computational study of glassy dynamics. 
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I. INTRODUCTION 



A deeper understanding of the behavior of disordered systems, and notably glassy ones, 
represents an outstanding challenge in modern condensed matter physics. In the face of the 
lack of a comprehensive theory, much work is currently devoted to the numerical simulation 
of glassy behavior. To date, these simulation efforts rely mainly on Molecular Dynamics 
and Monte Carlo techniques [1,2]. Both methods incorporate details of microscopic behavior 
(intermolecular potentials) and, as a result, they fall short of reaching spatio-temporal scales 
of macroscopic interest. A typical molecular dynamics simulation would last a few tens- 
hundreds of nanoseconds, to be contrasted with relaxation times of the order of an hour and 
more of real glasses. Larger scales are made accessible by coarse-grained lattice gas models 
in which the details of true molecular interactions are replaced by suitable heuristics about 
the behavior of fictitious lattice molecules. These models reach out larger scales but the 
corresponding outcomes have to be critically weighted against the simplifying assumptions 
they are based upon. At the extreme of this line one finds spin-glasses [3,4] and their modern 
lattice glass variants [5,6]. 

Although to a different degree of physical fidelity, all of these models retain the many- 
body nature of the intermolecular interactions, with an inevitable price in terms of com- 
putational effort. It is therefore of great interest to explore whether mesoscopic models, 
namely suitable effective one-body formulations may capture at least some basic features of 
glassy behavior, at a fraction of the computational cost associated with the aforementioned 
many-body simulation techniques. 

In the present paper this task is pursued within the framework of the lattice Boltzmann 
equation (LBE), a minimal kinetic equation describing stylized pseudo-molecules evolving 
in a regular lattice according to a simple local dynamics including free-streaming, collisional 
relaxation and (effective) intermolecular interactions [7-9] . In particular, we shall show that 
it is possible to construct a class of suitable mesoscopic potentials supporting some basic 
traits of disordered systems. Here we shall focus on a minimal requirement: The ability to 
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sustain sharp density contrasts within a disordered pattern. 

The paper is organized as follows. In the next Section, after a cursory view of the main 
features of the standard LBE, we present our model. Section III is devoted to the results 
of numerical simulations and finally, we discuss the main findings of the present paper and 
outline possible future research directions. 

II. THE MODEL 

The most popular form of lattice Boltzmann equation (Lattice BGK, for Bhatnagar, 
Gross, Krook) [10], reads as follows: 

fiif + c t ,t+ 1) - f t (f,t) = -u\Si - m(r,t) (1) 

where fi(r,t) = f(r,v = q,£) is a discrete distribution function of particles moving along 
with discrete speed q. The right-hand side represents the relaxation to a local equilibrium 
ff in a time lapse of the order of u; -1 . This local equilibrium is usually taken in the form of 
a quadratic expansion of a local Maxwellian: 

,e r, . U-Cj UU • (CjCj - cgl) 

fi = P^l 1 + — + ] ( 2 ) 

where p = Y,i fi is the fluid density, pu = Y^ifA the fluid current, Wi is a set of weights 
normalized to unity, c s the lattice sound speed (a constant equal to 1/ \/3 in our case) and 
finally, / stands for the unit tensor. The set of discrete speeds must be chosen in such a way 
as to guarantee mass, momentum and energy conservation, as well as rotational invariance. 
Only a limited subclass of lattices qualifies. In the sequel, we shall refer to the nine-speed 
lattice consisting of zero-speed, speed one c = 1 (nearest-neighbor connection), and speed 
c = a/2, (next-nearest-neighbor connection). Standard analysis shows that the LBE sys- 
tem behaves like a quasi-incompressible fluid with an ideal equation of state p = pc 2 s and 
kinematic viscosity v — c 2 s (1/uj — 1/2). Departures from ideal-gas behavior within the LBE 
formalism are typically described by means of phenomenological pseudo-potentials mimick- 
ing the effects of potential energy interactions [11]. These pseudo-potentials have proved 
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capable of describing a wide range of complex fluid behaviors, such as phase-separation 
in binary fluids and phase-transitions. To date, however, no LBE model appears to have 
addressed frustrated systems. This is precisely the scope of this work. 

We choose the empirical pseudo-potential in the form of a lattice-truncated repulsive 
Lennard- Jones interaction: 

V(r) = G^/(p)(r/a)- 12 , 1< r < V2 

V(r) = 0, elsewhere (3) 

where the cut-off a = 1 in lattice units and G is a parameter controlling the strength of the 
interaction (an inverse effective temperature), hence the surface tension between dense and 
light phases. The density-dependence is taken in the form of a hierarchy of polynomials, 
whose density zeroes are distributed according to a binary tree of depth N g (number of 
generations): 

Ng _ 1 29- 1 / _ \ 

*m=(p-po) n (1+r - WF -,) 2 n^f| (4) 

9=2 I 1 + ) k=l yP9k ~ PO) 

where 5po is the typical scale for the density gap. Starting with the root value po, at 
each generation g, a mother density p g k spawns two children Pg+i^k, Pg+i,2k+i spaced by an 
amount ±Sp /2 9 ~ 1 away from the mother value. This hierarchical potential implements the 
presence of N p = 2 Ng — 1 competing density extrema. 
As an example, for N g = 3, 

9 = 1 : Po 

9 = 2: pi, 2 = po(lTW2) (5) 
9 = 3 : p 3 ,4,5,6 = Po(l T Spo/2 =f 5p /A) 

The denominator serves the purpose of letting ^ — > outside the hierarchical range of zeroes 
p g k- The result is an effective potential alternating stable repulsive regions > 0, prime 
meaning derivative with respect to p) with unstable attractive ones < 0), distributed 
approximately in the density range p T = po(l =F Sp (l — 2~ Ns )). Since the case N g — 2 



(cubic interaction) is somewhat equivalent to a Van der Waals interaction, we can think 
to our potential as of a series of hierarchically nested van der Waals interactions. It is 
intuitively clear that if the coupling strength G is sufficiently strong, attractive regions are 
subject to instabilities which can be likened to a phase transition. However, since there are 
multiple forbidden regions interweaved with density gaps with (meta) stable behavior, the 
system is amenable to a full cascade of phase-transitions, depending on the strength of the 
interaction. As a result, the fluid is expected to exhibit a sort of "frustration" resulting 
from the competition of the multiple density minima. This competition is resolved through 
steep interfaces connecting the stable density regimes. In the language of free- volume theory, 
this frustration relates to the disordered short-range coexistence of two competing species: 
"voids" (regions with p < p ) and "cages" (regions with p > p ). This contrasts with binary 
separation scenarios, in which the two species organize over coherent patterns ("blobs") of 
sizeable extent. 

The LBE with hierarchical interactions takes the final form 

fi{f + 4 t+l)- fi(f, t) = -u[f t - /f] (f, t) + Fi[p(f, t)} (6) 

where Fj = — G-^V^[p] ■ Ci. This can be regarded as a dynamic, mean-field model of a 
non-ideal fluid with multiple density phases. 

III. NUMERICAL SIMULATIONS 

To assess the viability of the present LBE, we have simulated a 2D dimensional fluid 
with the following parameters: Grid size 256 x 256, u = 1.0, p = 1.0, N g = 3, 5p = 0.5, 
and G is varied between and 1.0. The resulting effective potential \I>(p) is shown in Fig. 1. 
The initial condition is p(x,y) — po (1 + £) where £ is a random perturbation uniformly 
distributed in the range [—0.01 : 0.01]. Due to the combined effect of kinetic hopping from 
site to site (the left-hand-side of the LBE) and mode-mode coupling [12] induced by the 
non-linear potential ^(p), the density distribution spreads out in time, so that an instability 
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is triggered as soon as G exceeds the critical threshold G c . In Fig. 2 we show the maximum 
and minimum densities as a function of the coupling strength G. This Figure highlights a 
sizeable symmetry breaking of the order of Ap/ p > 1 for G > G c , where Ap = p max — p m in- 
The density gap Ap starts up at G = G c ~ 0.09 and rapidly fills the available range of states 
between p_ and p + , as the interaction strength is increased. The scaling law is approximately 
Ap ~ (G — Gc) 1 / 2 in accordance with mean-field theory. The critical value estimated on the 
basis of the stability condition 

1*1 > cl/G (7) 
yields G c ~ 0.089, in a nice agreement with numerical data. A similar message is sent by the 



order parameter M = y< (5p) 2 >/ < p >, where 5p = p— < p > is the density fluctuation 
around its spatial average, which is represented by the dotted line in Fig. 2. Inspection of 
the spatial distribution of the density field p(x, y) shows that this non-zero order parameter 
is not associated with phase-separation, but rather with a disordered coexistence of voids 
and cages (see Fig. 3). This shows that the potential allows to model sharp density contrasts 
over a disordered spatial distribution. 

As a further observable, in Fig. 4 we plot the density-density correlation: 

< 5p(x + r x ,y + r y )5p(x, y) > 
9[r) <p(x,y)p(x,y)> ' W 

where r 2 = r 2 + r 2 . The curve has been obtained upon averaging the quantities g(r) over ten 
independent runs with different initial configurations. The error bars show the maximum 
errors on the values g(r) at each value of r. This function exhibits oscillating-decay behavior. 

Being in a lattice, it is clear that data are only available at discrete positions 1, y/2, 2, The 

small error bars and the non-oscillating behavior of g(r) in the absence of potential, which 
is virtually zero as a consequence of the uniform and constant density distribution, suggest 
that the observed oscillations in the case G ^ are not a numerical artifact. However, only 
further simulations of larger systems can confirm this result. 

In order to look for the existence of slow relaxation modes, we computed the density- 
density time-correlation function h{r) (see Fig. 5) 
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= < 5p(x,y,t + T)5p(x,y,t) > t 
< p(x,y,t)p(x,y,t)> t 

taken at x = 128, y = 128. In the above, < ... > t denotes a time average. We have run 
ten independent simulations starting from different initial configurations, and we computed 
h(r) in each run. These individual quantities were first averaged over time intervals of 150 
time-steps in order to eliminate high frequencies, and then the final average over the ten 
runs was taken. These quantities are plotted in Fig. 5 where the error bars represent the 
maximum error on the values of h{r). We find persistent oscillations around zero, possibly 
reflecting an everlasting competition between the various density minima of the potential 
over the timescale of the simulation. The case G = shows an immediate decay to zero of 
the time-correlation function, reflecting the fast relaxation towards a uniform distribution 
when no potential is acting on the system. 



IV. CONCLUSIONS 

We have presented a new mesoscopic lattice Boltzmann model with a density-dependent 
hierarchy of attractive and repulsive interactions. The long-time dynamics of this model 
yields disordered patterns with sharp short-range density contrasts, as opposed to binary 
phase-separation scenarios. Such configurations do not relax towards a frozen state, but 
continuously change preserving a disordered spatial structure. The present results show 
that one distinguishing feature of disordered behavior, namely frustration, can be modeled 
within the framework of a suitable one-body effective kinetic theory. On the other hand, 
the same results also show that further crucial signatures of glassy behavior, such as slow 
time-relaxation, are not captured by the present hierarchical model. At this point, the main 
question is: what major ingredients of glassy dynamics are we missing in the present model? 

According to Palmer et al [13], a successfully theory of glassy relaxation should satisfy 
three basic requirements: (1) The theory must be based on dynamics, and not just statistics, 
because since glasses break ergodicity, equilibrium distributions are of scanty use. (2) The 
above dynamics must be constrained, because there is no reasonable hope to diagonalize 
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such highly non-linear systems into independent modes. And finally, (3) the theory should 
be hierarchical, so that slow modes can adiabatically enslave the fast ones, while the fast 
modes, in return, set constraints on the dynamics of the slow ones. 

Our hierarchical LB model does indeed generically meet with all of these three require- 
ments, but only incompletely. The model is dynamic, but nonetheless based on a single-time 
relaxation to local equilibrium. Our model is also non-linear and hierarchically constrained, 
but still within the class of effective one-body theories. Such effective one-body theories 
prove exceedingly successful for simple fluids, but it might be that the physics of glasses, 
and particularly dynamic inhomogeneity [14,15], is just too intimately related to many-body 
effects which are beyond the reach of effective mean field theories, no matter how ingenious. 
We believe that, independently of the computational pay-offs discussed earlier in this work, 
this latter question bears a significant theoretical interest on its own. Work in this direction 
is currently underway. 
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FIG. 2. The maximum (full line) and minimum (dashed line) densities and the order parameter 
M (dotted line) as a function of the coupling strength G. 
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FIG. 3. Snapshot of the density configuration at t = 200 for the case G = 0.45. Black/white 
colors code for high/low density. 
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FIG. 4. Density-density correlation function at t = 200 for the cases G = (*) and G = 0.45 
(•), averaged over ten independent runs. Error bars are shown for the case G = 0.45. 
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FIG. 5. Density-density time correlation function for the cases G = (dashed line) and 
G = 0.45 (full line), averaged over ten independent runs. Error bars are shown for the case 
G = 0.45. 
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